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Abstract. The internal shock model for gamma-ray 
bursts involves shocks taking place in a relativistic wind 
with a very inhomogeneous initial distribution of the 
Lorentz factor. We have developed a ID lagrangian 
hydrocode to follow the evolution of such a wind and 
the results we have obtained are compared to those of 
a simpler model presented in a recent paper (Daigne & 
Mochkovitch 1998) where all pressure waves are sup- 



pressed in the wind so that shells with different velocities 
only interact by direct collisions. The detailed hydrody- 
namical calculation essentially confirms the conclusion 
of the simple model: the main temporal and spectral 
properties of gamma-ray bursts can be reproduced by 
internal shocks in a relativistic wind. 
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1. Introduction 

Since the discovery of the optical counterpart of GRB 



970228 (van Paradijs et al. 1997) the accurate localiza- 
tions provided by the Beppo-SAX satellite have led to 
the detection of the optical afterglow for more than ten 
gamma-ray bursts (hereafter GRBs). The most spectacu- 
lar result of these observations is to have provided a direct 
proof of the cosmological origin of GRBs. The detection 
of absorption lines at z — 0.835 in the spectrum of GRB 



970508 (Metzger et al. |1997|) followed by other redshift 
determinations (between z = 0.43 and z = 3.41) con- 
firmed the indications which were already available from 
the BATSE data showing a GRB distribution perfectly 
isotropic but non homogeneous in distance (Fishman and 
Meegan 1995| and references therein) . 

The energy release of GRBs with known redshifts ex- 



tends from E~, 



2 10 51 || to E-y 



2 10 54 ^r erg 



The 
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solid angle Q in which the emission is beamed is quite un- 
certain. A small fl should reveal itself by a break after a 
few days in the afterglow light curve. A break is indeed 
observed in a few cases such as GRB 990510 (Harrison et 
al. 1999) where could be as small as 0.01. However, 
most afterglows do not show any break which means that 
fl is usually not very small ~ 0.1 ?). 
The source of cosmic GRBs must therefore be able to re- 
lease a huge energy in a very short time. Possible candi- 
dates include the coalescence of two neutron stars (Eichler 



et al. 1989; Paczyhski 1991), the disruption of the neu- 



tron s tar in a neutron star - bl ack ho le binary (Narayan 
et al 



19921 Mochkovitch et al. 



1993 ) or t he col lapse of 
a massive star (Woosley 1993 , Paczyhski 1998 ). In all 
these cases the resulting configuration is expected to be 
a stellar mass black hole surrounded by a thick disc. 
Since the power emitted by GRBs is orders of magni- 
tude larger than the Eddington limit it cannot be radi- 
ated by a static photosphere. The released energy gen- 
erates a fireball which then leads to the formation of a 
wind. Moreover, this wind has to become highly relativis- 
tic in order to avoid the compactness problem and pro- 
duce gamma-rays (Baring 1995 ; Sari & Piran 1997 ). Val- 
ues of the Lorentz factor as high as T = 100-1000 are 
required, which limits the allowed amount of baryonic pol- 
lution to a remarkably low level. Only a few mechanisms 
have been proposed to produce a wind under such severe 
constraints : (i) magnetically driven outflow originating 
from the disc or powered by the Blandford-Znajek (1977) 

Mcszaros & Rees 1997; Daigne 
1999) ; (ii) reconnection 



process (Thomson 1994 
& Mochkovitch 



1999 



Lee et al. 

of magnetic field lines in the disc corona (Narayan et al. 
1992) ; (Hi) neutrino-antineutrino annihilation in a fun- 
nel along the rotation axis of the system (Meszaros & Rees 



1992; Mochkovitch et al. 1993, 1995). Mechanisms (i) and 



(ii) require that the magnetic field in the disc reaches very 
high values B > 10 15 G. Our preliminary study (Daigne 



& Mochkovitch 1999 ) of the wind emitted from the disc 
shows that it can avoid baryonic pollution only if some 
very severe constraints on the dissipation in the disc and 
the field geometry are satisfied. Some recent works (Ruf- 
fert et al. 1997) have also shown that mechanism (Hi) 
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is probably not efficient enough to power a gamma-ray 
burst, except may be for the shortest events. 

When the wind has reached its terminal Lorentz fac- 
tor, the energy is mainly stored in kinetic form and has to 
be converted back into gamma-rays. Two main ideas have 
been proposed to realize this conversion. The first one is 
the so-called external shock model (Rees & Meszaros 1992; 



Meszaros & Rees 1993). The wind is decelerated by the ex- 
ternal medium, leading to a shock. Gamma-rays are emit- 
ted by the accelerated electrons in the shocked material 
through the synchrotron and/or inverse Compton mecha- 
nisms. This model has been studied in details (Fenimore 
et al. 1997; Panaitescu et al. 1997: Panaitescu & Meszaros 



1998) and seems unable to reproduce some important fea- 



tures of GRBs such as their stro ng tem poral variability 
(see however Dermer & Mitman 1999 ). Conversely, the 
external shock model reproduces very well the delayed 
emissio n at l ower energy from the afterglows (Meszaros 
& Rees |l997t Wijers et al. |1997D . 

The sec ond p roposal is the internal shock model (Rees 
& Meszaros 1994 ) where the wind is supposed to be formed 
initially with a very inhomogeneous distribution of the 
Lorentz factor. Rapid parts of the wind then catch up 
with slower ones leading to internal shocks where gamma- 
rays are again produced by synchrotron or inverse Comp- 
ton radiation. We have started a study o f this model in 
a previous paper (Daigne & Mochkovitch 1998 , hereafter 
DM98) where the wind was simply made of a collection of 
"solid" shells interacting by direct collisions only (all pres- 
sure waves were suppressed). The very encouraging results 
we obtained had to be confirmed by a more detailed study. 
We have therefore developed a relativistic hydrocode to 
follow the evolution of the wind. We present the code and 
the main results in this paper. We write in Sect. 2 the 
lagrangian equations of hydrodynamics in special relativ- 
ity. In Sect. 3 we describe the numerical method we use 
to solve them and we present the tests we performed to 
validate the method. We display our results in Sect. 4 and 
Sect. 5 is the conclusion. 



2. Lagrangian equations of hydrodynamics in 
special relativity 

We write in a fixed frame the equations of mass, momen- 
tum and energy conservation in spherical symmetry : 



dV m _ d (R 2 v) 

at 

dS m 







dm 
am 



at 

dE m d (R 2 Pv 



at 



di 



= 



(1) 
(2) 
(3) 



where R and t are respectively the spatial and temporal 
coordinates in the fixed frame. The following quantities 



appear in Eqs. |^-|^ : P is the pressure in the fluid local rest 
frame, v is the fluid velocity in the fixed frame and V m , 
S m and E m are the specific volume, momentum density 
and energy density (including mass energy) in the fixed 
frame. These three quantities are related to quantities in 
the fluid local rest frame : 

(4) 



V, 



hTv , 



E m = hr-~ , 
r p 

where T — 



(5) 
(6) 

is the Lorentz factor, p is the rest-mass 



density and h = l + e + — is the specific enthalpy density (e 
being the specific internal energy density). The lagrangian 
mass coordinate m is defined by 



R- 



dR , 



(7) 



'-Rmin *™ 

-Rmin being the radius of the back edge of the wind. The 
system of Eqs [TJ-|] is completed by the equation of state 

P={ 1 -l)pe, (8) 

the adiabatic index 7 being a constant. 

In the non-relativistic limit (v — > and h — > 1), the 
quantities V m , S m and E m become equal to their newto- 
nian counterparts -, v and E = 1 + e + ^- and Eqs. [l]-|| 
then reduce to the classical equations of lagrangian hy- 
drodynamics in spherical symmetry. The great similarity 
between the relativistic and classical equations will allow 
us to use the powerful numerical methods which have been 
developed in classical hydrodynamics to follow the evolu- 
tion of a fluid with shocks. 

3. Numerical method 

3.1. Extension of the PPM to ID lagrangian relativistic 
hydrodynamics in spherical symmetry 

An extension of the Piccewise Parabolic Method of Colclla 



and Woodward ( 1984 ) to ID eulerian relativistic hydro- 
dynamics in planar sy mmet ry has been already presented 
by Marti and Muller ( 1996 ). We follow exactly the same 
procedure to extend the PPM to the lagrangian case in 
spherical symmetry. 

We adopt W = (V m , S m , E m ) as the set of variables. If 
the mass-averaged values WJ 1 of W at time t n in each cell 

m _i,m ,1 extending from R n _i to R n ,i arc known, 

J 2 J' 2 _ 3 2 J~r 2 

the values W™ +1 at time t n+1 = t n + At are computed in 
four steps: 

a) Reconstruction step 

arc obtained from W in each 



The variables U=(±,P,v 
cell by solving the following equation in h 



h 2 + (~f-l)h- 1 E m y/Sl + h 2 + jSi = 



(9) 
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Once h is known, the other quantities are easily computed 
since v = , Sm and - = rV m . These mass-averaged 

values are then interpolated by polynomials in the way 
descri bed in the original paper by Colella and Woodward 
( 1984 ). We use the same modifications of the coefficients of 
the interpolation polynomials, leading to a steeper repre- 
sentation of discontinuities and a monotone representation 
of smoother parts. 
b) Effective states 

At each interface j + h, two effective states U 4 



and 



Uj + i r (L and R denote the left and right sides of the 
interface) are constructed by mass-averaging these quan- 
tities in the region of cells j and j + 1 connected to the 
interface j + \ by a characteristic line during the time 
step. 

c) Riemann solver 

At each interface the two effectives states Uj + i L and 
Uj + i r define a Riemann problem (two constant states 
separated by a discontinuity surface), which is known to 
give rise, like in newtonian hydrodynamics, to two new 



states U. 



3+: 



,L* and U j+1 



separated by a contact dis- 
continuity and related to the initial states Uj + i L and 
Uj + 1 fj either by a shock or a rarefaction wave. The com- 
mon values of the pressure and the velocity of the two new 
intermediate states are given by the implicit equation 



= t'L*(P*) = "fl*(P*) ■ 

An analytic expression (for a polytropic gas) of 



vs*(p) 



IZsip) if P < Ps (rarefaction wave) 
Ssip) if p > Ps (shock wave) 



(10) 



(11) 



(where S either refers to th e L or R state) has been worked 
out by Marti and Miiller (|l994| ). Equation (ff) is solved 
using Brent's method to obtain the pressure and the 
velocity of the intermediate states at each interface. 
d) Time advancement 

The quantities Wf are calculated with numerical fluxes 
at each interface obtained from p* and v* in the same way 
than Colella and Woodward ( |f984| ). 

3.2. Numerical tests 

3.2.1. Relativistic shock tube 

We have sucessfully checked our code against two usual 
tests in relativistic hydrodynamics. The first one is the 
shock tube problem which is simply a Riemann problem 
in which the initial states are at rest. We present in Fig.[l] 
the results for pl — I, p L — 1000, vl = and pr = 1, 
Pn — 0.1, vr = 0. The adiabatic index is 7 = I and the 
discontinuity is initially located at x = 0.5. The figure is 
plotted at a time t = 0.303 for a grid of 1000 zones ini- 
tially equally spaced. The agreement between the exact 
and numerical profiles is satisfactory. The positions of the 
contact discontinuity and the shock are very accurate. The 
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Fig. 1. Relativistic shock tube problem with pi, — 1, pl = 
1000 and pu = 1, pu = 0.1 : Exact (solid line) and numer- 
ical profiles of density, pressure and velocity at t — 0.303. 

density, pressure and velocity of the post-shock state are 
also exact. However, the value of the density in the imme- 
diate vicinity of the contact discontinuity shows a small 
non-physical increase, which is more pronounced when the 
shock is stronger and disappears when the shock is weak 
(as in the shock tube problem with pi, = 10, pr, = 13.3 
and pr — 1, pr — 0, which has been considered by several 
authors). In the context of the internal shock model for 
GRBs, the shocks are only mildly relativistic and we do 
not observe any unexpected increase of the density in the 
results presented below. 

3.2.2. Spherical shock heating 

This test consists in a cold fluid, which is initially homo- 
geneous (p(R,0) = po) and enters a sphere of radius 1 
at constant velocity t>o- The fluid bounces at R = and 
is heated up. We present in Fig.|] the results for po = 1, 
Po = 10~ 6 and v Q = -0.99999 at t = 1.90. The adia- 
batic index is 7 = | and we used a grid of 1000 zones 
initially equally spaced. In the considered case of a cold 
homogeneous fluid, an analytical solution is known. The 



shocked state is at rest with a density p' = 

*• n 



7r +i 

7-1 



and a pressure p' Q = (7 — l)(r — 1)/5q- At time t the 
unshocked cold fluid of velocity vq has a distribution of 

density p(R,t) = po ^1 + ^r^J ■ The numerical profiles 

appear accurate except in the vicinity of the origin. The 
shock propagates with the correct velocity and the post- 
shock values of density and pressure are well reproduced. 
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Fig. 2. Spherical shock heating : Exact (solid line) and 
numerical profiles of density, pressure and velocity at t = 
1.90. 



We therefore conclude that the treatment of the geomet- 
rical terms in Eqs is correct. We have not tried to 
improve the computation near the center, which is not 
of major importance in the context of the internal shock 
model for GRBs where most of the emission takes place 
far from the origin. 

4. Results and discussion 

We have used our code to follow the evolution of a rel- 
ativistic wind with a very inhomogeneous initial distri- 
bution of the Lorentz factor. The first results have been 
already presented for small values of the Lorentz factor 
r ~ 40 (Daigne & Mochkovitch 1997). Here we describe 



our results for the large Lorentz factors (T > 100) which 
are relevant for the study of GRBs. We first consider the 
case of a simple single-pulse burst. 

4--1- Initial state 

Whatever the initial event leading to a GRB may be (NS- 
NS or NS-BH merger, "hypernova", etc), the system at 
the end of this preliminary stage is probably made of a 
stellar mass black hole surrounded by a thick disc (the 
"debris" torus). We consider that E, a substantial frac- 
tion of the available energy of the system, is injected at 
a typical radius i?o into a wind emitted during a dura- 
tion tw with a mass flow M. We do not discuss here the 
physical processes controlling M, tw and E but we as- 
sume that the baryonic load - — Mt ^ c is very small. 



Fig. 3. Initial state at t = 10 s for tw = 10 s. An energy 
E = 2 10 52 /4-7r erg/sr has been injected into the wind, 
whose mass is 7.76 1 28 g (which corresponds to an av- 
erage Lorentz factor T ~ 290J. The masses of the fast 
(T = 400 j and "slow" parts (T = 100 -> 400 J are equal. 
Upper panel : eulerian distribution of the Lorentz factor 
in the wind. Lower panel : corresponding lagrangian dis- 
tribution. 



The wind converts its internal energy into kinetic energy 
during its free expansion in the vacuum (the effect of the 
interstellar medium is negligible at this early stage) and 
accelerates until it reaches a Lorentz factor L ~ r/ at a 



typical radius TRq (Meszaros et al. 1993). This is where 
our simulation starts. 

More precisely, we define our initial state as follows. 
We consider that from t = to t = tw, a wind with a 
distribution of the Lorentz factor defined by 



_ I 250 150 cos (tt^ J if t < 0.4 t w 
\ 400 if t > 0.4 t w 



(12) 



has been produced by the source and that its back edge 
has reached i? m in = 400 Rq = 1.2 10 4 km (we adopt 
i?o = 30 km). We suppose that energy is injected at a 
constant rate E (we adopt E = 2 *° erg.s _1 /sr in the 
following), so that the total energy injected into the wind 
simply equals E = E tw and the injected mass flux is 
M(t) = jtjM c > 2 . The density profile in this initial state can 
be calculated if we assume that the internal energy is very 
small compared to the kinetic energy, which is indeed the 
case when the wind has reached its terminal Lorentz fac- 
tor (^p- <C 1). The eulerian and lagrangian profiles of L 
in the wind at t — tw are shown in FigJ| for tw — 10 s 
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Fig. 4. Left panel : Paths of the back and front edges and of the forward and reverse shocks in the t-t a plane. The 
two shocks appear at ts, the forward shock reaches the front edge very soon at 4f and the reverse shock reaches the 
back edge later at tfj. Right panel : Distribution of density p and Lorentz factor T at different times: (a) t = 2.5 10 4 s, 
just before the formation of the two shocks; (b) t = 3.0 10 4 s: the two shocks are clearly visible; (c) t = 3.1 10 4 s : the 
forward shock has just reached the front edge. The reverse shock has still more than one half of the mass to sweep; 
(d) t = 5.6 10 5 s: just before the reverse shock reaches the back edge. 



and E 



erg/sr. We have adopted 



1(T 3 and 



have checked that the results do not depend on this small 
value. 



4-2. Dynamical evolution 

The fast part of the wind catches up with the slower 
one. The matter is strongly compressed is the collision 
region, the velocity gradient becomes very steep and at 
ts ~ 3 10 3 tw, two shocks appear: a forward shock reach- 
ing the front edge at tp ~ 4 lO 3 ^ and a reverse shock 
reaching the back edge at tn ~ 9 I0 4 t\y- The hot and 
dense matter behind these two internal shocks radiates 
and produces the observed burst. The radiation losses are 
not taken into account in the dynamics, which is proba- 
bly not a too severe approximation since the dissipated 
energy represents about 10% of the total kinetic energy of 
the wind. 

When the two shocks have reached the edges, the evo- 
lution becomes unimportant regarding the emission of 
gamma-rays: two rarefaction waves develop at each edge 
and the wind continues to expand and cool. In fact, at 
this stage, the interstellar medium should absolutely be 
included in the calculation. An external shock propagates 
into the ISM which produces the afterglow and a reverse 
shock crosses the wind which can also lead to an observ- 



able emission. All these effects are not included in the 
present simulation, which is stopped when the two inter- 
nal shocks have reached the wind edges. 

We present in Fig.[| (left panel) the paths in a t - t a 
plot [t a = t — — is the arrival time of photons emitted 
at time t e = t on the line of sight at a distance R from 
the source) of the two shocks and of the two edges of the 
wind. In the right panel the corresponding distributions 
of r and p are plotted at different times. 

4-3. Gamma-ray emission and properties of the observed 
burst 

4.3.1. Method of calculation 

Consider an internal shock located at a distance R = R e 
from the source at a time t = t e in the fixed frame. The 
density p*,s, the Lorentz factor T^s and the specific in- 
ternal energy e*,s of the shocked (*) and unshocked (S) 
material are known from our hydrodynamical simulation. 
This shock will produce a contribution to the GRB which 
will be observed at an arrival time 

t a = t e ~ — (13) 



and which will last 

a Re 



2cTl 



(14) 
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Fig. 5. Single pulse burst (10s ) : Emission time i e , Lorentz 
factor of the emitting material IV, dissipated energy per 
proton tdiss and density of the shocked material p as a 
function of arrival time t a . Both contributions of the for- 
ward and reverse shocks are represented (the contribution 
of the forward shock is hardly visible). 

where T r is the Lorentz factor of the emitting material for 
which we adopt T r = T*. The luminosity of the shock is 
estimated by 



L s h — M s h (e» — es) 



(15) 



where M s h is the mass flux across the shock and tdiss = 
e* — eg is the dissipated energy per unit mass in the frame 
of the shocked material. 

Our code detects all the internal shocks present in the 
wind at a given time and saves their parameters in order 
to sum all the contributions to the emission and produce a 
synthetic gamma-ray burst. In a recent paper (DM98) we 
presented a simple model where the wind was idealized by 
a collection of "solid" shells interacting by direct collision 
only (i.e. all pressure waves were neglected). We detailed in 
this previous paper our assumptions to treat the emission 
of a given shock. We adopt here the same assumptions. 
The magnetic field in the shocked material is supposed to 
reach equipartition values 



B, 



(16) 



with oib = \- The Lorentz factor of the accelerated elec- 
trons is calculated using the expression given by Bykov & 
Meszaros (1996) who consider the scattering of electrons 
by turbulent magnetic field fluctuations : 

V(3-/*) 

(17) 
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Fig. 6. Single pulse burst (10s) : Magnetic field B eq , 
Lorentz factor of the accelerated electrons L e , synchrotron 
energy e syn and fraction of the energy which is radiated by 
the synchrotron process f syn as a function of arrival time 
t a . Both contributions of the forward and reverse shocks 
are represented (the contribution of the forward shock is 
hardly visible). 

where a m = 0.1 - 1 is the fraction of the dissipated 
energy which goes into the magnetic fluctuations; £ is 
the fraction of the electrons which are accelerated and 
A* (1-5 < M ^ 2) is the index of the fluctuation spectrum. 
For C ~ 1 an d M = 2, Eq. (17) corresponds to the usual 
equipartition assumption, leading to L e of a few hundreds. 
In this case, the emission of gamma-rays could result from 
inverse Compton scattering on synchotron photons. Bykov 
and Meszaros however suggests that only a small fraction 
C ~ 10~ 3 of the electrons may be accelerated, leading to T e 
values of several thousands. In this last case, synchrotron 
radiation can directly produce gamma-rays of typical en- 
ergy 



E syn = 500 



B, 



eq 



300 1000 G 



10 4 



keV 



(18) 



We now present a detailed comparison of the results of 
our hydrodynamical code with those previously obtained 
with the simple model (DM98) for a single pulse burst. 
We have plotted in Fig. || the values of t e , L, e diss and p 
as function of t a for the forward and reverse shocks. We 
observe an overall similarity between the two calculations, 
despite the crude approximations of the simple model. Not 
surprisingly, the worst estimated quantities are the post- 
shock density and the dissipated energy per proton, which 
are underestimated by a factor of ~ 5. Conversely, the 
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Fig. 7. Burst profiles for the initial distribution of the 
Lorentz factor shown in Fig. ^. The photon flux (normal- 
ized to the maximum count rate) is given in the interval 
50 - 300 keV, corresponding to BATSE bands 2+3. (a) 
Profile obtained with the expression of T e given by Eq. 
(|l7|); (b) same as (a) in logarithmic scale, which illustrates 
the exponential decay after maximum; (c) profile obtained 
with a constant T e ; (d) same as (c) in logarithmic scale. In 
the four panels, the full line represents the profile obtained 
with the hydrocode while the dashed line corresponds to 
the simple model. 

emission time and the Lorentz factor of the emitting ma- 
terial are correctly reproduced. The emission starts earlier 
in the simple model where there is no preliminary phase 
of compression before the formation of shocks (this leads 
to a larger underestimate of the density at the very begin- 
ning of the simulation), and ends later. The total efficiency 
of the dissipation process is also smaller ~ 5% instead of 
12% for the detailed model. 

The other quantities B eq , T e and e syn are not di- 
rectly given by the hydrodynamical simulation but arc 
parametrized by as, aw, £ and /i, whose values are un- 
known. To make a useful comparison between the two se- 
ries of results, we take the same as and fi in the two 
cases but adjust ckm/C 80 that the typical synchrotron en- 
ergy is the same. The corresponding values of B eq , T e and 
e syn are represented in Fig. ^| with as = 1/3, /i = 1.75 
and «m/( = 100 for the hydrocode and 1000 for the sim- 
ple model. As expected because of the differences in den- 
sity and dissipated energy, the magnetic field is under- 
estimated by a factor of 5 in the simple model. This is 
corrected by our choice of parameters for T e and the re- 
sulting synchrotron energies are very similar in the two 



Fig. 8. Ratio of the decay to rise times as a function of 
burst duration. The initial distribution of the Lorentz fac- 
tor is given by Eq. ( |l2| ) and the total energy injected into 

the wind is proportional to tw ■ E = 2 |S (to~s) er g/ sr - 
The full line corresponds to the results of the hydrocode 
while the dashed line shows the same relation obtained 
with the simple model. 

cases. Also notice that the efficiency of the synchrotron 
process is smaller in the simple model due to a poor esti- 
mate of the mass flux accross the shock. 

The agreement between the two calculations is satis- 
factory and allows to be quite confident in the results of 
the simple model. Compared to the hydrodynamical code, 
the simple model has very short computing times and en- 
ables a detailed exploration of the temporal and spectral 
properties of synthetic bursts which was presented in our 
previous paper (DM98). We show in the next section the 
detailed results obtained with the hydrocode in the case 
of a single pulse burst. 

4.3.2. Temporal properties 

The contributions of the forward and the reverse shocks 
are added to construct the synthetic burst. We assume 
that the photons emitted from t to t + dt by an internal 
shock of current luminosity are distributed according 
to a simple power-law spectrum 

d(En(E)) L sh dt ( E \~ x 

— If — x ~F ~F~ ' ( 19 ) 

where we adopt x = 2/3ora; = 3/2 (the two extreme low 
energy index that are expected for a synchrotron spec- 
trum) for E < E syn and 2 < x < 3 for E > E syn (x = 2.5 
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Fig. 9. Evolution of the profiles with duration. The initial 
distribution of the Lorentz factor is given by Eq. ( p"2| ) and 
the total injected energy is E — 2 *° {jo^) erg/sr. The 
different profiles correspond to (a) tw = 10 s ; (b) tw = 5 s 
; (c) tw — 2 s ; (d) tw = 1 s ; (e) tw — 0.5 s ; (f) tw = 0.2 s 
; (g) t w = 0.1 s. 



in the following). We take into account cosmological ef- 
fects (time dilation and redshift) assuming that the burst 
is located at z = 0.5. 

We have plotted in Figj^a the photon flux observed in 
BATSE bands 2+3 for the initial distribution of Lorentz 
factor shown in Fig. 3, calculated either with the hy- 



drocode (with = 100) or the simple model (with ° 



1000). The two profiles look similar but the hydrodynam- 
ical code leads to a slower decay. With i 5 (resp. £95) being 
the time when 5% (resp 95 %) of the total fluence has been 
received, we obtain a duration T90 = £95 — ts = 10.4 s in- 
stead of 6.67 s with the simple model. Figure illustrates 
that the exponential decay of the burst is also nicely repro- 
duced with the detailed calculation. However, if we define 



tmax as the time of maximum count rate and t t = t r , 



-U 



and Td = t< 



95 



t, 



as the rise and the decay times, we 
get a ratio Td/r r = 2.08. DM98 found that a larger value 
of Td/r r and a corresponding profile closer to the charac- 
teristic "fast rise - exponential decay" (FRED) shape is 
obtained by assuming that the fraction £ of accelerated 
electrons increases with the dissipated energy per proton 
ediss- As in DM98 we adopt £ oc edi SS , so that T e is inde- 
pendent of ediss ■ Figures 0c and f?]d show the resulting pro- 
files with T e = 5000 for the hydrodocode and T e = 10000 
for the simple model. The profile then better reproduces 
a typical FRED shape. 

The observed ten dency of short bursts to become sym- 
metric (Norris et al. 1996| ) has been tested in DM98 with 




100 

E (keV) 



1000 



Fig. 10. Spectrum of the burst presented in figure |^(c) 
(tw = 2 s). The number of photons per energy inter- 
val n(E) and the product E 2 n(E) are shown in ar- 
bitrary unit. This product is maximum at peak energy 
E p = 403 keV in case (a) (x = -2/3) and E p = 193 keV 
in case (b) (x — —3/2). The dashed lines show a fit of each 
spectrum with Band's formula in the interval 10 keV - 10 
MeV (parameters are given in the text). 



the simple model. The basic behaviour was reproduced 
but the effect was even exagerated since, for T go < 1 s, 
i~d/T r was smaller than unity i.e. the decline was faster 
than the rise. As can be seen in Fig. || and ^, the situation 
is improved with the hydrocode since now Td/r r ~ 1 for 
T 90 ~ 0.4 s. However, the shortest bursts are still asym- 
metric with Td/r r ~ 0.6 for T90 < 0.2 s. 

Figure || also shows that the ratio Td/r r is limited to 
a maximum value of ~ 2.5 for the longest bursts which 
appears to be in contradiction with the short rise times 
observed in some cases. As discussed in DM98, an initial 
distribution of the Lorentz factor with a steeper gradient 
than the one used here (Eq. (12)) can indeed increase Td/r r 
but extreme values (such as Td/r r possibly larger than 10 
in GRB 970208) might still be difficult to reproduce. 

4.3.3. Spectral properties 

In DM98 we presented a complete study of the global and 
instantaneous spectral properties of synthetic bursts cal- 
culated with the simple model. Since these spectral prop- 
erties are hardly different when calculated with the hy- 
drocode, we do not present them in detail again. We just 
show in Fig. |Io| the shape of the global spectrum calcu- 
lated for the single pulse burst. Despite the very simple 
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form adopted for the instantaneous spectrum (Eq. p^ ), 
the sum of all the elementary contributions produces an 
overall spectrum with a more complex shape, which is well 



reproduced with Band's formula (Band et al. 1993) 



n(E) 
n{E) 



A 



A 



E 



100 keV 

(a - 0)E O 



r.xj> J ~ ) for E •_ (n - ,')£•„ 



Q-/3 



100 keV 

for # > (a - 0)E O 



exp (a — /?) 



100 keV 



(20) 



We find values of the parameters comparable to those ob- 
served in real bursts. The best fits in Fig. [l^ correspond 
to a = -0.935, = -2.42 and E = 239 keV in case (a) 
(x = -2/3) and a = -1.60, = -2.47 and E = 609 keV 
in case (b) (x = —3/2). As x is limited to the range 
2/3 < x < 3/2, we cannot get spectra with low energy 
slopes flatter than —2/3 as they are observed in several 
bursts (Preece et al. 2000 ). A more detailed description of 
the radiative processes is then needed to reproduce these 
extreme slopes (an attemp t to s olve this problem is pro- 
posed by Meszaros & Rees 2000 ). 

However, even with the crude modelization of the instan- 
taneous spectra which is used here, it has been shown in 
DM98 that several spectral properties of GRBs are re- 
produced. In particular, the hard to soft evolution during 
a pulse and the change of pulse shape as a function of 
energy as well as the duration - hardness ratio relation 
which appears as a natural consequence of the internal 
shock model. These important spectral features are con- 
firmed in our detailed hydrodynamical calculation. 

4-4- Case of more complex bursts 

An important property of the internal shock model is its 
ability to produce a great variety of temporal profiles. Nor- 
ris et al. ( 1996[ ) have shown that complex bursts can gen- 
erally be analysed in terms of a series of (possibly over- 
lapping) simple pulses. This result is readily interpreted in 
the context of the internal shock model. A wind made of a 
succession of fast and slow shells will produce a succession 
of pulses which will add to form a complex burst. 

We present such examples of complex bursts in Fig. 11 
and 12. The first one (Fig. [ll]) is produced by an initial 
distribution of the Lorentz factor made of five consecu- 
tive identical patterns. Each pattern made of a slow and 
a rapid part produces its own individual pulse and the re- 
sulting burst has a complex shape with five, very similar, 
pulses. Our second example (Fig. 12) uses the same type 
of initial distribution of the Lorentz factor but the slow 
parts now have non equal T values. The resulting burst is 
more realistic with four pulses of different intensities. 

We did not treat with the hydrocode a large number 
of cases as we did with the simple model. Nevertheless 
we confirmed the essential result that the variability in- 
troduced in the initial distribution of the Lorentz factor 



is present in the burst profile with the same time scale. 
The profile therefore appears as a direct indicator of the 
activity of the central engine. 



5. Conclusions 

This paper is the continuation of our study of the inter- 
nal shock model started in DM98. We developed a ID la- 
grangian relativistic hydrocode (in spherical symmetry) to 
validate our previous simpler approach where all pressure 
waves were neglected in the wind. Our code is an extension 
of the classical PPM method of Colella and Woodwards 



( |1984D in the spirit of the work by Marti & Miiller ( |1996D 
for the eulerian case in planar symmetry. 

A detailed comparison has been made between the 
hydrocode and the simple model in the case of a single 
pulse burst. It appears that the dynamical evolution of 
the wind is well reproduced by the simple model, which is 
not too surprising because the wind energy is largely dom- 
inated by the kinetic part so that the effect of pressure 
waves is small. Only one physical quantity - the density 
of the shocked material - is strongly underestimated in 
the simple model. In order to make valuable comparisons 
between the two calculations we have therefore adjusted 
the equipartition parameters so that the mean value of 
the synchrotron energy is the same in the two cases. The 
synthetic bursts which are then obtained are very simi- 
lar which proves that our first approach was essentially 
correct and confirm our previous results. A similar con- 



clusion was reached by Panaitescu and Meszaros (1999) 
who performed a comparable study. 

The internal shock model can easily explain the great 
temporal variability observed in GRBs. The main charac- 
teristic features of individual pulses are well reproduced: 
(1) pulses have typical asymmetric "FRED" profiles; (2) 
the pulse width decreases with energy following a power- 
law W(E) cx E~ v with p ~ 0.4; (3) short pulses show 
a tendency to become more symmetric. Our model still 
gives very short pulses which decay faster than they rise 
but the hydrodynamical simulation improves the situa- 
tion compared to the simple model. Spectral properties of 
GRBs are also well reproduced. We obtain synthetic spec- 
tra which can be nicely fitted with Band's function with 
parameters comparable to those observed in real GRBs. 
The spectral hardness and the count rate are correlated 
during the evolution of a burst with the hardness usually 
preceeding the count rate. As also pointed in DM98, the 
duration-hardness relation is a natural consequence of the 
internal shock model. These results are very encouraging 
and the main difficulty which remains is the low efficiency 
(about 10%) of the internal shock model. As long as the 
energetics of GRBs and the mechanism initially operating 
in the central engine are not precisely identified, we can- 
not say if this is a critical problem or not. We still believe 
that the internal shock model is at present the most con- 
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Fig. 11. Example of a complex burst. Left panel: initial distribution of the Lorentz factor with five identical patterns. 
Right panel: corresponding profile obtained with the hydrocode (normalized photon flux in the interval 50 - 300 keV). 
The wind is produced with a duration tw — 10 s and the injected energy is E = 10 52 /47r erg/sr. 




Fig. 12. Another complex burst. Left panel: initial distribution of the Lorentz factor with now four non identical pat- 
terns (Lorentz factors in the slow parts are different). Right panel: corresponding profile obtained with the hydrocode 
(normalized photon flux in the interval 50 - 300 keV). The wind is produced with a duration tw — 10 s and the 
injected energy is E = 10 52 /47r erg/sr. 



vicing candidate to explain the gamma-ray emission from 
GRBs. 

Next steps in this work will address the following ques- 
tions. We first want to extend our hydrodynamical code to 
a non-adiabatic version in order to include the radiative 
losses in the dynamical calculation. We have already de- 
veloped an "isothermal Rieman Solver" for that purpose 



(Daigne & Mochkovitch 1997| ) . We would also like to study 
the effects of the external medium, with a special atten- 
tion to the reverse shock which propagates into the wind 
and possibly interacts with the internal shocks. Prelimi- 
nary results with the simple method using "solid layers" 



have already been obtained (Daigne & Mochkovitch 1999) 
but they have to be confirmed by a hydrodynamical calcu- 
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lation. Finally, we would like to investigate the details of 
the emission process during internal shocks to solve some 
of the problems encountered by the synchrotron model. 
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